Spatial fluctuations of a surviving particle in the trapping reaction 
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»' ' We consider the trapping reaction, A + B B, where A and B particles have a diffusive dynamics 

characterized by diffusion constants Da and Db- The interaction with B particles can be formally 
{2J[J' incorporated in an effective dynamics for one A particle as was recently shown by Bray et al. 

I [Phys. Rev. E 67, 060102 (2003)]. We use this method to compute, in space dimension d = 1, 

the asymptotic behaviour of the spatial fluctuation, {z^{t))^^^, for a surviving A particle in the 
perturbative regime, Da/Db 1, for the case of an initially uniform distribution of B particles. 
We show that, for t ^ 1, {z^(t))^^^ oc t''' with (j> = 1/4. By contrast, the fluctuations of paths 
constrained to return to their starting point at time t grow with the larger exponent 1/3. Numerical 
tests are consistent with these predictions. 
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I. INTRODUCTION 



Reaction-diffusion problems are a paradigm of irreversible nonequilibrium dynamics 0. Among the conceptually 
simplest of such problems is the trapping reaction, A + B B, in which A and B particles diffuse in space. When 
, an ^ and B particle meet the A particle is annihilated, while there are no interactions between particles of the same 
Ch • species. Though simple to state, this problem is quite subtle: for example, the asymptotic behaviour of the A-particle 
I [ density for one- and two-dimensional systems has oiily recently been computed [3, 13 . 

. Recently there has been renewed interest 0, IE IE H IOj HJ in this reaction, stimulated in part by the advent 
C ■ of both new numerical methods and new analytical techniques 0, Q . The main focus of interest in this system 
^ , has traditionally been the computation of the time-dependence of the ^-particle density ^E|- Note that, since the 
A-particles are independent, this computation is equivalent to the computation of the survival probability, Q{t), of a 
single ^-particle moving in a sea of B particles. In the case of one space dimension, to which this paper is devoted, 
Q{t) is known to decay asymptotically as exp(— Ai^/^) Recently, by a combination of bounding arguments, 

^ ] the constant A was determined exactly |EIE|- = (2/7r)/9(47rZ?B)^/^, where p and Db are the density and diffusion 
constant respectively of the _B-particles. Note that the constant A does not depend on the A-particle diffusion constant, 
Jyi-^ Da, and the leading asymptotic result is the same as for a stationary A-particle. 

■ A parallel advance in algorithmic methods has enabled very long simulations of the model, especially in one 
f*^ ■ dimension 0,113, in which Q{t) can be made arbitrarily small (of order 10 , say, in reasonable computer time). 
Remarkably, the predicted asymptotic decay is not observed in simulations even for these very small values of Q{t), 
motivating a study of the preasymptotic behaviour. In a companion paper |l4j . we have employed a path-integral 
approach j^j to compute, for small ratio Da/Db of diffusion constants, the leading correction to the asymptotic form 
quoted above. Including this term gives reasonable agreement with data. 

In the present paper we address, within the path-integral formalism, a different aspect of the problem. Specifically, 
. we study the spatial fluctuations, {x'^{t)), of the yl-particle averaged over the surviving trajectories at time t. We find 
^ ' that the fluctuations are subdiffusive: {x^{t))^/^ ^ with tf) < 1/2. In their numerical work, Mehra and Grassberger 
O ■ 03 i^n*i 4> to be in the range 0.25 — 0.3, where the effective exponent seems to be decreasing at larger t. We will argue 
O' that 0=1/4. 

Our approach is based on a recent development in ref. ■ The main idea is to consider the probability of intersection 
before time t, P{x,t; z), between a B particle starting from x at t = and a fixed trajectory, z{t), of the A particle 
starting at x = when t — 0. This probability can be averaged over the initial position of the B particle (assumed 
uniform in length V, where we specialize to space dimension d = 1): R{t;z)/V = (1/V^) / dxP{x,t; z). The survival 
probability, i.e. the probability of no collision up to time t, for the trajectory z of the A particle is 1 — R{t] z)/V. If 
the system has N _B-particles, uniformly distributed in the length V , the survival probability for the given trajectory 
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is (1 — R{t; z)/V)'^ which, in the continuum Umit N — > oo, V oo, with p = N/V = constant, becomes exp{— 
where p,[z] — pR(t; z). 

The statistical properties of the A particle are calculated by averaging over the trajectories, where each trajectory is 
weighted with its survival probability multiplied by its Brownian weight. In this way one obtains the action functional 

i 

which defines an effective dynamics for surviving A particle trajectories 0|. 

Further progress requires an explicit expression for the functional p\z\. In ref. :9i| it was shown that the constant 
trajectory, z{t) — 0, gives the global minimum of and furthermore a systematic expansion in Da/Db was 

obtained for ii[z\. To first order one obtains |9j 

= At^/^ + ^ /* 7r4W r 1/2, ,,,J ^(^i) - -(^^r + 0{z') , (2) 
8t:Db Jo [t-hV'^ Jo t^\ti-t2Y/^ 

where A = (2/7r)p(47rDB)^/^- 

In this paper we study the asymptotic properties of the spatial fluctuation of the A particle, using the formalism 
described above in the quadratic approximation for given by Eq. The paper is laid out as follows. In 

section II we apply the formalism developed in refs. 9, 14] to compute the root-mean-square fluctuation of surviving 
trajectories, and show that it scales as i^/^. In section III we adapt this formalism to compute the root- mean-square 
fluctuation at time t/2 for surviving 'loop' trajectories, constrained to return to their starting points at time t. This 
second calculation is motivated by the fact that it involves a very different mathematical treatment, and it is not 
obvious that the same fluctuation exponent will be obtained. Indeed, we find that this second quantity scales with a 
larger exponent 1/3, which seems at first sight counter to physical intuition. Numerical results presented in section III 
are, however, consistent with these two different predictions. Section IV is a summary and conclusion. Some technical 
details are relegated to an appendix. 



II. SPATIAL FLUCTUATIONS OF SURVIVING TRAJECTORIES 



The probability for the A particle, starting from the origin and following a trajectory z{t)^ to reach x at time t can 
be written as 

p{x, t)^N r^'^ " V[z{t)] exp(-^[z]) . (3) 

Since we are interested in the spatial fluctuations of the surviving trajectories we choose the normalization factor M 
such that 

dxp{x, t) ^ 1 . (4) 

We can extract the x-dependence of the functional integral, Eq. Q, by expanding the trajectories z{t) around the 
'classical path', Zdit), which minimizes the action S. Since the action is a quadratic functional of z{t), we have 

l.zit)=0 

p(x,t)^NeyiY>{-S[za-x,t]) exp(-5[z]) , (5) 

J z{Q)=0 

and the remaining functional integral, over the fluctuations around the classical path, can be absorbed into the 
normalization constant, as can the z-independent part of see Eq. 

To compute the probability distribution p{x, t), therefore, it is sufficient to flnd the classical path Zc/, i.e. the path 
of least action. The analysis is simplifled by the introduction of a dimensionless action, S, through the substitutions 
ti tu, t2 tv, z{t) — > xS, in Eq. jSJ. This gives 

s^'-^s^u' Miuf + Ut) r r - av)r , (6) 
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where 7 = {\/2T:){DA/DB)t^''^ = {-n"^ /2)g{t) is a dimensionless constant, the final equality defining g. In what follows, 
it is important to bear in mind that g oc t^/'^ , and that we are interested in the large-t behaviour of the theory. 

The minimum action can be computed using the observation that any trajectory starting at ^(0) = and ending at 
^(1) — 1 can be written as a straight line connecting the end points plus a Fourier expansion for the deviation from 
this straight line, 

^(u) — u + a„ sin(n7rM) , (7) 

ri>l 

where the linear part u can also be written as a Fourier sine series: 

^^^_2^sin(n^^ (8) 

n>l 

With these substitutions we can express the action, Eq. JHJ, in terms of the matrix elements Amn introduced in |l4j . 

1 du r dv 



mn 



y/1 -u Jq y/v{u~ v)^/'^ 



X [sin(n7r7i) — sm(m:v)] [sin(m7r?i) — sin(m7rti)]. (9) 
After a straightforward computation we obtain: 

1 TT^ 

" 9 + S^Anri - 1^9 ^ A^nnan + — ^{n^5mn + .9-4„i„nm)a,„a„. (10) 



2 -y / J — I itti ,j / J — iiLit — ft ' ^ 

m,n m,n m,n 



The classical path is obtained by minimizing this expression with respect to the coefficients a„, which gives 

TT ^-^ \ I + qA / nm 

m,g 

where we have made the rescaling a„ = a„/n. Hence the minimum action can be expressed in terms of the matrix 
elements Amn- 



mn rn,n ^ 



pq-^mn 



= ^ + 5E(7^)„.>- (12) 

The desired probability distribution p{x, t) has the Gaussian form 

(X exp(-S'd) = cxp ^-^^^x^^ , (13) 

with variance (x^) = DAt/Sd- We emphasize at this point that Sd depends on time through g. Our next goal, 
therefore, is to determine the g-dependence, at large g, of the minimum action given by Eq. (|12|l . To this end we 
employ a finite-size scaling method, by computing numerically the action Sd for an A matrix truncated to size N, 
and analyzing the dependence of Sd on g and N. 

To motivate the particular finite-size scaling form we use, we first consider the extreme limits g — > 00 and g — > 
at fixed N. For g 00, Eq. ifT^ becomes Sd = 1/2 -I- J2n=i I = N + 1/2, while for g ^ we obtain Sd = 
1/2 + gJ2n m=i^""i + O(g^). From the definition ^ of A„in it is easy to show that the double sum converges 
to the value for N 00. Thus we can write Sd — 1/2 = H{g,N), where H{g,N) — > N for g — > 00 and 

H{g, 00) {TT^/Q)g for small g. This suggests the finite-size scaling form Sd — 1/2 ~ Nh{g/N), where h{x) ~ x for 
small x so that the A^-dependence drops out as iV ^ 00. 

A scaling plot, testing this form, is presented in Figure 1 for a range of values of g and A^. The data collapse 
is excellent (the data sets are virtually indistinguishable) confirming the validity of the assumed scaling form. We 
conclude that, in the limit N 00, Sd ^ g ^ t^^^, implying, via Eq. (|13|l . that 

(zHt)) o, p-\DBty/^ . (14) 
Note that {z'^{t)) is independent of Da at leading order. 
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g/N 

FIG. 1: Scaling of the g-dependent part of the dimensionless extremal action for free trajectories, for N = 50, 100, 200, 400. 
The continuous straight line with slope 1 is a guide to the eye. It indicates the expected asymptotic slope of the data for 
g/N < 1. 



III. FLUCTUATIONS OF CLOSED LOOPS 



We turn now to a discussion of the spatial fluctuations of closed loops, i.e. trajectories that return to the origin at 
time t. One might naively expect that the root-mcan-square displacement of the trajectory at, say, its mid-point t/2, 
would also scale as t^^'^. The mathematics of this calculation is, however, rather different and we shall see that the 
final result is different as well. 

The probability distribution of the displacement x of such a trajectory at time t/2 is given by 



p^^'^'P^ {x, t/2) = M / V[z{t)]S{x - z(V2))e-^W . 



(15) 



Applying the standard 'exponentiation' of the delta function via an auxiliary integration, representing z by a Fourier 
sine series on the interval (0, t), and performing the resulting Gaussian integrals we obtain, after some straightforward 
algebra, 



DC exp 



2DAt 



Fig) 



(16) 



where 



St 



9 A 



\{-Km/2) sin(7i"n/2) 



mn 



(17) 



It follows that 



{x\t/2)) 



DAt _ 4:0 At 



n,n=l ^ ^ 



sin(7rrn/2) sin(7rn/2) 



mn 



(18) 



Note the very different character of the calculations of p{x,t) and p^^°°P\x,t/2). The former is completely deter- 
mined by the classical path connecting the spacetime points (0,0) and {x,t). The fluctuations around this path are 
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FIG. 2: Scaling of the function F{g,N) (see text) for closed trajectories at time t/2, for A'' — 50, 100, 200, 400: continuous 
curves - full problem; broken curves - diagonal approximation. The straight line with slope 2/3 is a guide to the eye. It 
represents the expected asymptotic slope of the data for g/N^''^ <^ 1. 



not important for the calculation of p{x, t) since they give an x-independent contribution, i.e. they only the affect 
the normalization of p{x,t), which can in any case be fixed by hand a posteriori. By contrast, the calculation of 
p^^°°p)[x^t/2)^ over paths connecting the spacetime points (0,0) and (0,i), is entirely concerned with fluctuations - 
we need to explicitly integrate out the fluctuations to determine p'^^°°P\x , t /2) . 

To determine {x'^ {t/2)) loops from Eq. (|16|l we need to determine the form of the function F{g) for large g. Again 
we can exploit a numerical finite-size scaling analysis. In view of the considerations in the preceding paragraph, we 
expect a rather different finite-size scaling form for the function F(g,N), obtained from truncating the matrix A to 
size N, than we found for S{g, N). In this context it is helpful to recall the calculation of the subleading term in 
the A-particle survival probability, which is also a fluctuation-dominated effect. In that calculation it was found that 
approximating the A matrix by its diagonal part gave a very good guide to the finite-size scaling behaviour. Using 
that approach here, we are led to consider the function 



Fd^agig,N) = 




For large n, Ann n '^/^ (see the appendix) and, for large N and g it is clear that large n dominate the sum, justifying 
the use of the large n form of Ann- Inserting this form in Eq. H19|) we see that, for TV oo, F^iag ~ g^^^ esc t^^^ for 
large g, whence Eq. ((TB|l gives {x^ {t/2)) loops ~ F)At/g'^^^ cx p~^^^{DADB)^^^t'^^^ within the diagonal approximation. 

Furthermore, the diagonal approximation H19|) exhibits, for g and N both large, the scaling form Fi^agig, N) — 
g'^/^H{g/N^I'^) = Nh{g/N^^^). We adopt this same scaling form for the analysis of the full function F{g,N), 
determined numerically. The results, together with the diagonal approximation, are shown in Figure 2. The data 
collapse is not as good as for the full set of paths (Fig. 1), with the data drifting to higher values for larger N. The 
same is true, however, of the diagonal data for which we know by explicit calculation that the assumed asymptotic 
scaling form is correct. Indeed, the full function F{g,N) appears to be converging faster with increasing N than the 
diagonal approximation to it. We are therefore confident that the assumed scaling is asymptotically correct and, as 
a consequence, a/ {x"^ {t/2)) loops t^^^. 

An intriguing aspect of this result is that a subset of A-particle trajectories have spatial fluctuations scaling with a 
larger exponent {{x"^)]/^^ ~ t^^^) than the whole set of trajectories ((cc^)^/^ ~ t^^'^). While this result seems at first 
sight counterintuitive, the numerical scaling analysis strongly suggests that it is correct. Further supporting data will 
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be presented in the following section. 

It is worth noting that, in this context, there is nothing special about the fluctuations of the mid-point of a closed 
path. For example, we can consider the mean square fluctuation, (a;^);^^^^ — t^^ dt' {x^{t')) loops, along the path. 
Using the usual Fourier expansion one readily obtains 

n=l ^ ' y / nn 

This quantity will also scale as t^^^ (i.e. the root-mean square fluctuation grows as t^^^) as can be checked explicitly 
within the diagonal approximation. As a further twist, we can consider the fluctuations at time t/2, or the root-mean 
square fluctuation over the whole path, of all trajectories that survive until time t. These both scale at t^^^ in contrast 
to the fluctuation of the end-point of the path which, as we have seen, grows as 

That the fluctuations of the trajectory at time t/2, for paths which survive to time t, could be larger than the 
fluctuations at time t/2 of paths which survive till time t/2 has the following intuitive explanation. Trajectories which 
survive for a long time arise predominantly from initial configurations of the B particles in which the A particle is 
initially in a large region devoid of B particles. If the A particle is to survive until time i, this region is larger than 
if it only has to survive until time t/2, so the fluctuations in its position at time t/2 are expected to be larger in the 
former case than in the latter. 



IV. NUMERICAL SIMULATIONS 



Numerical simulation data were obtained for both the unconstrained average {z^{t))^/^ and the average {z^ {t / 2)) loops 
corresponding to surviving trajectories that return to the origin at time t. The data were generated using the 
lattice-based algorithm of Mehra and Grassberger which is readily adapted to the constrained case. The 

basic algorithm calculates the survival probability, averaged over all i?-particle initial conditions and trajectories, 
of a given A-particle trajectory. Final results are obtained by averaging over a large number (typically 10^) of A- 
particle trajectories randomly generated from the 2* possible paths of length t. In the modified algorithm, A-particle 
trajectories are generated dynamically. At a given time, the next step is taken to the left with a probability equal to 
the number of paths that reach the origin at time t ('return paths'), given that the first step is to the left, divided 
by the total number of return paths from that point. This procedure generates the ensemble of closed paths with the 
required weights. 

The parameter values used were Da = Db = 1/2 and p = 1/2. In both cases we divide the root-mean-square value 
of the displacement by t^^'^ and by t^^^ to test the different predictions for the two types of average. The data are 
presented in Figures 3 and 4. For ease of comparison, we have defined the end-time of the unconstrained paths to be 
t/2, i.e. all data are taken at the same real time. 

Looking at Figure 3, we see that the closed-path (loop) data show an increasing trend with time t consistent with an 
exponent larger that 1/4. A linear fit on a log-log plot gives an exponent ~ 0.33, consistent with the predicted value 
of 1/3. The unconstrained data also show an increasing trend, but there is a suggestion that the slope is decreasing 
at later times - the prediction = 1/4 requires that these data eventually reach a plateau as i ^ oo. A linear fit, 
on a log-log plot, gives ~ 0.30 over this range. The companion Figure 4, in which the data are scaled by t-^^^, 
tells a similar story. The closed-path data are consistent with (p = 1/3, but the data including all paths show a clear 
decrease with time, indicative of a smaller exponent. 

To explore this behaviour in more detail, we have computed a time-dependent effective exponent, (/>(t). For the 
computation we use the full data set, of which only a small subset (of the 'all paths' data) is displayed in Figs. 3 and 
4. The total time interval is divided into consecutive bins of size dti. The average value of {z^{t))^/'^ is computed 
in each bin to smooth the data prior to taking a numerical derivative. The effective exponent, (fiiti), at time ti 
is approximated by the difference between the logarithms of the average value for bins i and i + 1 divided by the 
difference of the logarithms of the average positions of bins i and i -t- 1. The result of this procedure is shown in 
Fig. 5. The effective exponent is noisy, but is clearly decreasing with increasing time in a manner consistent with the 
predicted limiting value of 1/4. The equivalent data for the closed paths is also shown. In this case, however, there 
are no more data available beyond that shown in Figures 3 and 4 - the simulations are computationally expensive due 
to the requirement that the paths return to the origin at time t. As a result, each bin contains a single data point, 
and the effective exponent is necessarily very noisy, but consistent with = 1/3. 
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FIG. 3: Root mean square displacement at time t/2, divided by t^^'*, for all paths and for closed paths that survive till time t. 
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FIG. 4: Root mean square displacement at time t/2, divided by t^^^, for all paths and for closed paths that survive till time t. 



V. CONCLUSION 



We have studied the scaling properties of the spatial fluctuations of a diffusing particle interacting with diffusing 
traps, using the effective action arising from eliminating the trap degrees of freedom. The presence of the trapping 
process renders the fluctuations in the position of surviving particles subdiffusive, (2;^(t))^/^ oc t^/*, with an amplitude 
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1/t 

FIG. 5: Effective exponent (f>{t) for all paths and for closed paths, plotted against 1/t. 



that is independent of the j4-particlc diffusion constant, Da- We have also found that the fluctuations at 'half-time' 



of surviving trajectories that return to the origin at time t scale with a different exponent {z'^{t/2)) 



1/2 
{loops) 



OC 



As 



a cautionary note we recall that both of these results were obtained by using the quadratic approximation to the 
action functional, which is strictly valid only in the limit Da Db- It remains an open question to what extent 
corrections to the quadratic approximation can modify the results for general values of Da/ Db- In view of the 
reasonable qualitative agreement of the theory with numerical data for the case Da = Db, however, we expect any 
modifications to be quantitative rather than qualitative. 
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APPENDIX A: THE MATRIX A 



The matrix element Amn is proportional to the double integral 



mnAr, 



dx 



dy 



[sin(n7rx) — sin(7?,7ry)] [sin(m7ra;) — sin(m7ry)] 



(Al) 



which, by the variable substitutions z = x — y, s — y and standard manipulations of the trigonometric functions can 
be reduced to a sum of one-dimensional integrals, 



27rcos 



2tt cos 



7r(TO — n)\ dz fmTZ\ / 



7r(m + \ dz 



o(|(m-n)(l-z)) 



f dz / Tin Z \ / TTITTZ \ / TT \ 



(A2) 



where Jo(^) is the Bessel function of the first kind. The fact that = 0ifm-|-n = 2p+l can be obtained directly 
from Eq. IjAljl on noticing that the kernel is symmetric on reflection about the line y ~ 1 — x. 
For the case m = n, the dominant large-rt contribution comes from the first integral: 



2n 



^ dz .2 /"mrz , _ 
_ sm ) « 2ttV^ 



dz 
372 



sm 



(A3) 



More elaborate is the calculation of the asymptotic behaviour for Imn at large n with m fixed and large. We do 
the same substitutions, z = x — y, s = y,m order to isolate the most singular part in Eq. (jAl|) . Then we integrate 
once by parts with respect to z such that the singularity in z becomes integrable: z^'^/^ z^^/^. In the next step 
we transform the products of trigonometric function into sums to give 



1 



-{n — to) Jn,„+m ~ {n — m)Jn^m,n-7n + + 'm)Jn+m,n+m + {n + m)Jn^n- 



where the integrals Jn,m and Kn,m are given by 



-m,n— m 



+ (n + m)Kn 



+ (n - m)Kn,n- 



K 



^ dz 



ds 



Tj2 ^1/2(1 „^)i/2 sin7r(nz + m(l-z)5) 
^ dz ds(l-2s) . , ^ N 

.1/2(1 _,)l/2 ^^"^(-^ + "^(1-^)-^) 



(A4) 

(A5) 
(A6) 



Now we have to calculate these integrals to leading order in the n ^ 1, m ^ 1 and n ^ m. We show in detail 
the computation for Jn,n+m.- We deform the integral from (0,1) to two parallel lines of the form — > ±ioo and 
1 ± ioo 1, the sign being chosen such that the intermediate segment ±ioo — > 1 ± ioo vanishes. This gives 



(TV, 

n,n+rn ^ 



ds 



^ dz 



Z^^^ Jo sl/2(l_s)l/2 

^ dz 



^i7T{nz-\-{n-\-m){l — z)s) 



(A7) 







^1/2- 

p27r7r(n+m) (1 — 2) 



ids 



Visl/2(1 _ i5)l/2 

ids 



-7r{n-\-m){l — z)s 



\/=Isl/2(l+is)l/2 
1 ^^giTrnz gi7r(n+m) 



— 7r(n+m)(l — £)s 



^/(n + m) Jo ^^(1 - z) Vi\/ {n + m) Jo v^z(l - z) 



^i7T(n+m) 



^(n + m) \ Vi 



Vi^/{n + m) I \M 



•^n + TO V ^/n 



C0S7r(n + to) 



TO 



(A8) 

(A9) 

(AlO) 
(All) 
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FIG. 6: Aran OS s functioii of n at m = 50, 55 starting from the main diagonal of the matrix. In the inset Ai„ is compared 
with the asymptotic result. The symbols are the from the numerical integration of Eq. lAlL the lines are the asymptotic 
approximation, Eq. I|A13|I . 



where we have used, to leading order, 

when (5(0) exists and fc ^ 1. 

Special attention must be paid to the point z = 1 in the integral over s, as the argument in the exponential is not 
large when z 1 with to, n fixed. However, one can show that this region gives a subleading contribution. 

Applying the same procedure to the remaining integrals, we find that the contribution of Jn-m,n-m, Jn+m.n+m 
and Kn-m,n-nn -f^n+m,n+m cancel each other to leading order and we finally obtain the closed- form expression 

(A13) 

mn \ -^n — m y to + n / 

for TO 3> 1, S> 1 and n ^ m. It is worth noting that the above result is a good approximation for any value of 
TO, n outside the main diagonal as one can see from Fig. This fact is useful in numerical computations since the 
integrals that define Amn have strongly oscillating integrands for large m, n. 



